其他
Newton-Raphson Method估算函数的根
10年前写的，数值计算还是有必要练一练的
Newton-Raphson Method
曲线f(x)有根c,取曲线上一点(x_1,f(x_1)), 过此点的切线交x轴x_2,过曲线上(x_2,f(x_2))的切线交x轴x_3，如此反复得到一个序列 x_1,x_2, …,x_n 逼近c值.
过(x_n,f(x_n))的切线方程为 y-f(x_n) = f’(x_n),(x-x_n),假设此方程与x轴的交点为x_n+1, 即有: 0 - f(x_n) = f’(x_n)(x_n+1 - x_n)， 即x_n+1 = x_n - f(x_n)/f’(x_n) <Eq. 1>.
下面利用此法来求一个数的开方。 f(x) = x^2 - a 有根sqrt(a),
由f’(x_n) = 2x_n, 代入式<Eq. 1>可得x_n+1 = (x_n + a/x_n)/2; 当i -> INF 时, x_i -> sqrt(a);
C implementation
#include <stdio.h>
int main(void){
int a,error;
double x0,x1 = 1;
do {
printf("Input a positive integer: ");
scanf("%d",&a);
if (error = (a <=0))
printf("\nERROR: Do it again!\n\n");
}
while (error);
while (x0 != x1) {
x0 = x1; /* save the current value of x1 */
x1 = 0.5 * (x1 + a / x1); /* compute a new value of x1 */
}
printf("%lf\n",x1);
return 0;
}R implemantation
2010-01-11 用R来实现一下，不单是求开方，估算函数的根。
newton <- function(fun, x0, tol=1e-7, niter=100) {
fun.list = as.list(fun)
var <- names(fun.list[1])
fun.exp = fun.list[[2]]
diff.fun = D(fun.exp, var)
df = list(x=0, diff.fun)
df = as.function(df)
for (i in 1:niter) {
x = x0 - fun(x0)/df(x0)
if (abs(fun(x)) < tol)
return(x)
x0 = x
}
stop("exceeded allowed number of iterations")
}> f = function(x) x^2 – 5
> newton(f, 4)
[1] 2.236068
> g = function(x) x^3 – 5
> newton(g, 4)
[1] 1.709976往期精彩